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ABSTRACT 



f*"*- ' One of the largest uncertainties in understanding the effect of a background 



UV field on galaxy formation is the intensity and evolution of the radiation field 
with redshift. This work attempts to shed light on this issue by computing the 



quasi-hydrostatic equilibrium states of gas in spherically symmetric dark matter 
. halos (roughly corresponding to dwarf galaxies) as a function of the amplitude of the 

background UV field. We integrate the full equations of radiative transfer, heating, 
cooling and non-equilibrium chemistry for nine species: H, H + , H~,H2, H^, He, He + , 
He ++ , and e~ . As the amplitude of the UV background is decreased the gas in the core 
\Q ' of the dwarf goes through three stages characterized by the predominance of ionized 

(H + ), neutral (H) and molecular (H2) hydrogen. Characterizing the gas state of a 
t^J- ■ dwarf galaxy with the radiation field allows us to estimate its behavior for a variety of 

models of the background UV flux. Our results indicate that a typical radiation field 
, can easily delay the collapse of gas in halos corresponding to 1-er CDM perturbations 

\ with circular velocities less than 30 km s _1 . 

Or 

o : 

£3 . Subject headings: cosmolgy:theory — galaxies:formation — molecular processes — 

radiative transfer 
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1. Introduction 

How do galaxies form? Most models predict that galaxies are assembled through a successive 
series of mergers of smaller systems, a process known as hierarchical clustering. These models also 
predict that small galaxies form before big galaxies. Observations, however, suggest the opposite: 



while low mass galaxies are forming the bulk of their stars at z ~ 1 (Driver et al 1995; Babul & 



Ferguson 1996| ), large galaxies appear to be well-established by z ~ 1 and do not show significant 



evolution between z ~ 1 and today in luminosity, color, size or space density (pteidel & Dickinson 



1994 |Steidel, Dickinson fc Persson 1994j) . 
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The resolution of this apparent conflict between theory and observations may lie in the 
physics of galaxy formation: the hierarchical clustering paradigm only describes the formation 
of gravitationally bound entities, not the process of converting gas into stars. Observations of 
Lya clouds along quasar lines-of-sight suggest that at high redshifts, the universe is permeated 
by a metagalactic UV flux that may suppress star formation by rapidly dissociating atomic 
and molecular gas. At z = 2, this background is estimated to have a strength J u ~ 10~ 21 
s 



erg 



1 cm 2 ster 1 Hz 1 at the Lyman limit ( [Bechtold et al 1987 ). As discussed by Babul fc Rees 



1992| (see also Efstathiou 1992 ), an ionizing flux of such intensity can easily prevent the gas in 



low mass halos from settling and forming stars. Analytic calculations ( Rees 1986 ; Ikeuchi 
suggest that the halos affected will be those with circular velocities v c < 30 km s _1 , a result 
confirmed by recent numerical studies ( |Katz, Weinberg Sz Hernquist 19"9€ ; Navarro & Steinmetz 



1996] |Thoul fc Weinberg 1996|) . 



Although the epoch-dependence of the UV flux is uncertain, it is expected that J u declines 
with redshift. A declining J u has two effects on the gravitationally bound, photoionized gas in 
the halos. First, the rate at which the gas is photoheated will decrease and hence, its equilibrium 
temperature will also decrease slightly. The reduced temperature results in a gradual concentration 
of the gas towards the center ( Ikeuchi, Murakami fc Rees 1989| ). Second, the gas becomes more 
neutral and offers a greater optical depth to the ionizing photons leading to a diminution of the 
ionizing flux reaching the central regions and the formation of a warm (T ~ 9000 K) shielded 



neutral core ( Murakami fc Ikeuchi 1990 ). The formation of a neutral core, however, is not a 
sufficient condition for star formation. For a gas cloud to be susceptible to star formation, it must 
be at least marginally self-gravitating (e.g., Matthews 1972; ) and this can occur only if the gas in 
halos with v c < 30 km s _1 can cool below 9000 K. 

Further cooling of a neutral metal-poor gas requires the formation and survival of molecular 
hydrogen. The formation of molecular hydrogen in a gas of primordial composition occurs via gas 
phase reactions. Various studies have shown that an external photoionizing radiation field can 
greatly affect the efficiency of formation of H2 molecules (Shapiro & Kang 1987; Haiman, Rees 



fc Loeb 1996a| ; |E aiman, Rees &: Loeb 1996b|) , with even a moderate UV flux being capable of 
suppressing the H2 abundance. In time, however, the decline in the UV intensity will result in 
the formation of sufficient molecular hydrogen to allow the gas in the halos to cool and become 
susceptible to star formation. Babul & Rees (1992) argued for z = 1 as the epoch of galaxy 
formation in low mass halos, and sought to identify these forming galaxies with the numerous 
small faint blue galaxies seen in deep images. 

In this paper, we present a more thorough investigation of the epoch of galaxy formation 
in smaller halos with circular velocities in the range 15 km s _1 < v c < 70 km s . Specifically, 
we only consider halos that form after the metagalactic UV flux has been established {z < 5). 
Through a detailed study of the thermal and ionization structure — we take into account radiation 
transfer, as well as heating, cooling and the corresponding non-equilibrium chemistry (Abel et 
al 1996| ; [Anninos et al 1996| ) for nine species: H, H + , H~, H2, H^, He, He + , He ++ , and e~ 
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of the hydrostatic-equilibrium configuration of halo gas subject to a range of UV intensities, we 
determine the threshold UV intensity at which molecular hydrogen begins to form. Given a 
specific model for the evolution of the UV flux, the threshold intensity can be easily converted 
into a redshift of galaxy formation. 

In §2 of this paper, we describe the dwarf galaxy model, which includes: the dark matter 
halos, the amount of gas expected in these halos, the details of the radiative transfer, heating, 
cooling and non-equilibrium chemistry, and the numerical method. §3 Discusses the results of our 
simulations. §4 Gives our conclusions. 



2. Dwarf Galaxy Model 

The current numerical models can not yet incorporate the wide range of physical processes 
and physical scales associated with galaxy formation. We focus on the quasi-hydrostatic evolution 
of a gas cloud in a small dark matter halo. By simplifying the problem to the time evolution of a 
spherically symmetric cloud in a dark matter halo, we were able to incorporate the non-equilibrium 
chemistry and treat the radiative transfer in much greater detail. 



2.1. Dark Matter Halo 

We will simulate the evolution of the gas in a fixed halo potential. For our purposes, the dark 
matter halo is specified by two parameters: the circular velocity v c and the virialization redshift 
z v , which can be translated into a halo radius rhaio and halo mass Mh a io if we assume that the 
overdensity at virialization is S ~ 175 ( |Gunn Sz Gott 1972| ) 

2 GM halo 4vr 3 
v c = , irr hSLlo 5p c (z) = M ha i , (1) 

r halo " 

where the mean density is given by the usual expressions for a 0, = 1 CDM cosmological model: 

Pc (z) = (1 + Zfpl OTrG^Lbble = 1, bubble = 2/3# . 

We will focus on the evolution of "typical" halos, which correspond to 1-a perturbations. In a 
b = 1.6 standard CDM model, Press-Schechter theory implies that the halos virialized at a redshift 
given by 

1 + z v « (7.5/1.7)(M halo /10 9 M ) ai , 10 8 M < M halo < 10 10 M (2) 

Figure 1 shows the relationship between halo circular velocity and the virialization epoch for the 
1-a perturbations. The basic trend is generic to all hierarchical models: small objects form first 
and larger objects form later. 



We use a dark matter halo profile that has been fit by Burkert 1995 to galaxy rotation curves, 



PDM{r)= (l + x)(l + x*) ' X = r/r ° (3) 
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which in turn can be related to the halo radius and mass by 

nmio = 3.4r , 
-Mhalo = M DM (rh a io), 



M DM (r) = / p DM (r)^r 2 dr. (4) 



o 



While recent numerical work suggests that the halo density profiles of large galaxies are 
proportional to r~ l in the centers and r -3 at the edges ( |Navarro, Frenk &; White 1996| ), these 



profiles do not fit the dwarf galaxy observations ( Moore 1994 ; Flores fc Primack 1994| ). 



2.2. Gas Content and Profile 

The ultraviolet background is able to heat the gas to a temperature of roughly 10 4 °K. In 
large halos, where v c > 50 km s _1 , the gas pressure is relatively unimportant and the gas content 
is determined by the global value of f^: M gas = f^Mhaio (assuming ft = 1). However, for smaller 
halos collapsing out of a hot IGM the gas pressure resists the collapse ([Thoul fc Weinberg 199(f ) 
and M gas < ^&Af na i . We now make some simple estimates as to where this transition occurs 
and how much gas should reside in the halo. If the gas in the uncollapsed halo is greater than 
the Jeans mass than the gas should collapse of its own accord. This provides an upper limit to 
amount of gas in the halo 

M gas = n 6 M ha i , IWhalo > Mj 

cans (5) 

where M Jeans = /lV/io c (27rciGM^Hubbic) 3 , / = 3/^^, Cj 2 GM = l.hk B Ti GM / For 
OfcMhaio < -/Vfj eans we can formulate an upper limit for M gas by estimating the amount of 
mass that could be acreted via Bondi accretion in a Hubble time. Thus, 

M gas = M Bo ndi , fWhalo < Mj (6) 

where M Bon di = /iH u bbieM Bondi , M Bondi = vrG 2 M^ alo /) c /c 3 GM . Note, the 0(1) factor / has been 
included in M Jeans and M Bondi so that M Jeans = M Bon di when 0, b M halo = M Jeans . Plots of these 
various masses as function of redshift are shown in Figure 2 for Tjgm = 20, 000 °K, indicating 
that the halos become more DM dominated as they get smaller, which is consistent with the 
observations ( Carignan Freeman 1988| ; |de Blok McGaugh 1996 ). 



The initial gas density profile is specified by hydrostatic equilibrium and by our assumption 
that the gas is in thermal equilibrium with the background radiation field. This approximation 
will hold as long as the gas has time to react to changes in the DM potential, self-gravity, and 
cooling. The scales governing these processes are the sound speed, Hubble time, Jeans mass, and 
the net cooling/heating time. The sound speed of a 10 4 °K gas is ~ 15 km s _1 , corresponding 
to a crossing time for a typical r na i ~ 10 kpc halo of ~ 6 x 10 8 years. The DM potential will 
change on the order of a Hubble time, which at z ~ 3 is ~ 2 x 10 9 years. The gas will be stable 
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to self-gravity if it is less than the Jeans mass, which from Figure 2 is true for nearly the entire 
range of models. The heating/cooling time becomes important only when a significant amount of 
H2 has formed, and determining this point one of the goals of this paper. 



2.3. Non-Equilibrium Chemistry 



The important role of the detailed chemistry of primordial gas (in particular the formation of 
H2) has been known and studied since it was first proposed as a mechanism for the formation of 
globular clusters (Peebles & Dicke 1968). The potential number of reactions in this simple mixture 
of H and He is enormous ( |Janev 1987 ). Abel et al 1996| have selected a subset of these reactions 
to model the behavior of primordial gas for low densities (n < 10 4 ) over a range of temperatures 
(1°K < T < 10 8 °K); these equations (see Appendix A) represent a careful balance between 
computational efficiency and accuracy. 



2.4. Radiative Transfer 

Fully 3D radiative transfer requires estimating the contribution to the flux at every point 
from every other point along all paths for each wavelength. At the minimum this is a 6D problem. 
However, in most instances symmetries can be introduced which result in a more tractable 
situation. The simplest situation occurs when the gas can be assumed to be optically thin 
throughout. This approximation is sufficient in the majority of cosmological situations ( [Katz, 



Weinberg & Hernquist 1996; Navarro & Steinmetz 1996; Anninos et al 1996) and only breaks 



down in the cores of halos that have undergone sufficient cooling, a situation that is usually made 
further intractable by the complexities of star formation. The next simplest geometry is that of a 
slab (or a sphere under the assumption of a radially perpendicular radiation field), which leaves 
an intrinsically 2D problem. This approach is the most common in radiative transfer and has 
been used to address similar situations ( Haiman, Rees fc Loeb 1996a ; Haiman Loeb 1996j ), 



Although this approach may not be a bad approximation for a sphere in an isotropic radiation 
field, we choose to account for all the different paths that penetrate a given spherical shell, which 
leaves an inherently 3D problem. Although this geometry can significantly increase the size of the 
computations, for a static grid (i.e. non-adaptive) pre-computing of the various geometric factors 
can significant alleviate this situation (see Appendix B). Taking into account the different paths 
effectively "softens" the optical depth, smoothing out transitions from optically thin to optically 
thick regimes. In addition, introducing the accounting for the different paths lays the groundwork 
for exploring more general geometries. 
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2.5. Heating and Cooling 

Perhaps the most important aspect of the model is the balance between the heating and 
cooling processes. This balance is what allows the establishment of a quasi-static temperature 
profile for a specific radiative flux. If the balance between the heating and cooling is not 
established, then the hydrostatic equilibrium solution to the gas profile will evolve too rapidly. 
Fortunately, this situation only comes about when the gas in the halo becomes dense and a large 
amount H2 is formed. As was mentioned earlier, this point presumably marks the onset of star 
formation, which is one of the central goals of this work. 



al 1996 



The temperature profile is evolved via the heating and cooling functions found in Anninos et 

T = T(T) - E(T) 
T ET { ' 

where E = 1.5kBTJ2 n i, i = e_ > H, H + , H~, H2, H^, He, He + , He ++ . T includes photoionizaton 
heating and E includes cooling due to collisional excitation, collisional ionization, recombination, 
molecular hydrogen, bremsstrahlung and Compton cooling. All the appropriate functions are 
taken from Appendix B of Anninos et al 1996[ 



2.6. Numerical Method 

The microphysical processes couple to the larger scale density profile primarily through 
radiative heating, which sets the temperature profile. The rate of radiative heating is in turn 
strongly dependent on the column densities of each species, which is set by the temperature. Thus, 
we have a system which is described by differential equations on the small scale with integral 
constraints on the large scale. The difficulty of solving such a system in the large variety of time 
scales involved. Solving the entire set simultaneously is prohibitive. Our approach has been to use 
code modules which solve for each of the processes independently. Iterating between the modules 
then provides adequate approximation to the true solution. The species solver and heating/cooling 
modules were provided by Yu Zhang and Mike Norman of NCSA ( [Abel et al 1996 ; Anninos~et| 



al 1996j ). The radiative transfer and hydrostatic equilibrium modules we wrote ourselves. Our 



central goal is to see how the gas behaves as a function of the amplitude of the background UV 
field. The overall approach is to pick a halo, set an initial value of Jo, and estimate the time and 
number of steps to integrate at this value. We then use each module at each time step which 
converges on the overall solution. Each module was tested and verified independently. In addition, 
where possible, combinations of the modules have been tested. A more detailed description of our 
numerical method is given in Appendix C. 
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3. 



Results and Discussion 



The goal of the simulation is to determine the maximum ionizing flux that permitted the 
formation of molecular gas (and stars). The simulations began with the gas ionized. Over the 
range of objects, the state of the core exhibits the same qualitative behavior (see Figure 3). As 
the flux decreases the core goes through three phases ionized H + (T corc ~ 20, 000°K), neutral H 
(T corc ~ 10,000°K) and molecular H 2 (T corc ~ 100°K). 

In the initial H + state, the object is completely ionized and resembles the IGM. In the H state 
the core is neutral and the object is like an inverse Stromgren sphere, except that the ionization is 
primarily collisional, a result of the high temperature maintained by the radiative heating. In the 
H2 state, conditions allow for the formation of molecular hydrogen in the core; consequently, the 
cooling time becomes much less than the dynamical time and the core collapses. 

Five simulations were conducted for halos in the range 15 km s _1 < v c < 70 km s _1 
corresponding to: 3 x 10 8 Mq < M halo < 5 x 10 10 Mq, 4 > z v > 2 , and 
3 x 10 5 Mq < M gas < 3 x 10 9 Mq (see Figure 2). For each object there are fluxes Jh,h 2 
at which the upper transition H + — > H and the lower transition H — ► H2 occur (see Figure 4). 
Jh,k 2 ( z v) can be fit by a formula 



where a H ,H 2 = (17.9,17.3), and 6 h ,h 2 = (6.2,6.4) 

Knowing Jh,h 2 allows us to place bounds on the behavior of an object in a given radiation 
field. Most significantly, if the radiation field is above Jh, then it will be impossible for the object 
to collapse and form stars. If the radiation field is below Jh 2 , then it must collapse. Most likely 
the actual value of the flux at which collapse occurs is between Jh and Jh 2 . 

For objects consistent with l-cr perturbations in a standard CDM cosmology selecting the 
evolution of the radiation field (J(z, v)) specifies critical redshifts zh,h 2 such that for zn 2 < z < zu 
the object will be in the neutral H state. Figures 5, 6 and 7 plot the virial redshift and the critical 
redshifts as a function of the circular velocity for three different amplitudes of an evolving UV field 



These results indicate that a typical radiation field can easily prevent gas in halos with 
v c < 30 km s _1 from collapsing. Furthermore, the less massive halos, which are typically the 
first to virialize, are the last to form galaxies. Thus, if these results have any correspondence to 
observed dwarf galaxies, they suggest that the larger objects will burst earlier, and perhaps will 
show an observable trend of increasing age with mass. 



ln[JH,H 2 /<^2l] = OH,H 2 — ^H,H 2 ^', 



•V ) 



(8) 




z < 3 
z > 3 



(9) 



Another important aspect of these objects is their potential contribution to Quasar absorption 
lines. Figure 8 shows the column density evolution for H and He + . During the H + phase 
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A^h ~ 10 16 ~ 17 , suggesting that these objects may make a contribution to the Lya forest and Ly 
limit systems. During the H phase Ah increases dramatically to levels near those seen in damped 
Lya clouds, (N^ ~ 10 21 ) but this is probably a relatively short lived phase and therefore is it 
unlikely that these objects would be seen in Quasar spectra. 

Our calculations ignored the role of metals in the evolution of the gas. Metals can have 
a number of important effects on the chemistry and dynamics of gas clouds: dust grains will 
absorb ionizing radiation and serve as formation sites for molecular hydrogen and atomic lines 
of C and other heavy elements can be important coolants at temperatures around 10 4 °K. Are 
these processes important in dwarf galaxies at high redshift? Observations of QSO Lyman forest 
clouds suggest that the metal abundances in meta-galactic gas was Z ~ 0.001 — 0.01 times the 
solar value at z ~ 3 ( [Songaira fc Cowie 1996] ) . At these abundances, heavy element cooling is 
unimportant ( [Bohringer fc Hensler 1989| ). If we make the conservative assumption that most 
of the carbon at high redshift is incorporated into dust grains (as in our galaxy), with a size 
distribution that is similar to the local ISM, then this suggests a cross-section per hydrogen atom 
of (TdustUOOOJl) ~ Z 2 x 10~ 21 cm 2 ( Draine &: Bertoldi 1996] ). In our models, the maximum 



column density occurs when the cloud is most centrally condensed, which occurs just before the 
onset of H2 formation and is roug hly AT H ~ 10 21 cm" 2 (see Fi gure 8). Thus the maximum optical 
depth at these wavelengths is approximately Td U st ~ Z ~ 0.01. The contribution of dust to H2 
formation in our galaxy can be approximated by a nn 2 

« Rnn H , where R = 6 x lO^^T 1 / 2 cm s 



( Draine fc Bertoldi 1996 ). If we scale R by the metallicity, then the dust term will be negligible 
in comparison to the other terms contributing to H2 formation whenever njj- > 10~ 9 , which is 
nearly always the case in the neutral H core. 



4. Conclusions 

In this paper, we have studied the evolution of a proto-dwarf galaxy exposed to the 
metagalactic radiation field. We begin with an ionized gas cloud, initially Jeans stable, in a 
dark matter halo and follow its hydrostatic evolution as the meta-galactic flux decreases. Our 
calculations indicate that the state of the gas can be characterized by the predominance of ionized, 
neutral, or molecular hydrogen in the core. The transitions between these phases takes place 
quickly as the amplitude of the background flux decreases. We have computed the critical fluxes 
at which the transitions take place Jh,h 2 > which serves as upper and lower bounds on the flux at 
which rapid cooling and the subsequent star formation can begin. If the flux is greater than Jh, 
then is unlikely that the gas will cool. If the flux is less that Jh 2 , then the gas must cool. 

Characterizing the state of the gas in terms of the background flux allows us to use any model 
for the evolution of flux with redshift. Using typical values of the flux, our simulations indicate 
that gas in a v c = 30 km s _1 halo collapsing at z ~ 3 can easily be prevented from forming 
significant amounts of H2 until z ~ 2. Thus our calculations are consistent with photoionization 
delaying the formation of low mass galaxies flBabul fc Rees 1992j ). 
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In this paper we have focused on 1-a perturabtions, in subsequent work we will expand our 
exploration of the (v c ,z v ) plane, which will allow us to further address the issues pertaining to 
the observed faint blue galaxies. As an example, our preliminary calculations indicate that a 
v c = 30 km s _1 halo which virializes at z v = 2 will first be able to form neutral H at z ~ 1.25. 



We gratefully acknowledge Yu Zhang and Mike Norman for the use of their non-equilibrium 
species solver. J. Kepner was supported by the Dept. of Energy Computational Science Fellowship 
Program. DNS is supported by the MAP/MIDEX program. This work was partially supported 
by NASA through Hubble Archival Research Grant AR-06337.20-94A awarded by Space Telescope 
Science Institute, which is operated by the Association of Universities for Research in Astronomy, 
Inc., under NASA contract NAS5-26555. 



A. Non-equilibrium Chemistry 

The evolution of the electron, hydrogen and helium abundances are computed from the 
following set of coupled equations ( Abel et al 1996| ): 



n H = k 2 n e -n n + - {k\n e - + k 20 )n-a (Al) 

n H + = -htt (A2) 

"He = hn e -n He + - (k 3 n c - + k 21 )n He (A3) 

«He+ = -n H e-n He ++ (A4) 

™He++ = k 5 n e -n He - k 6 n e -n Ke ++ + k 22 n He + (A5) 

_ fc 7 n e -n H ^ A6 ^ 

(h + fci5)riH + (fcie + k n)n K + + k M n e - + k 23 

n H2 = k 8 n H -nu + kion^+nn + ki 9 n u +n R - 

- (hsm + hin-u+ + ki 2 n c - + k 2A + k 27 )n R2 (A7) 

fc 9 n H nH+ + fciin H2 n H + + /ci 7 n H -n H+ + (k 24 + k 27 )n n2 , . Q . 

n H + = : : : : : (A8) 

2 honu + h 8 n e - + fci 9 n H - + k 25 + k 2e 

h e - = (fcxnn - k 2 n n + + k 3 n Ue - /c 4 "- H e+ + hn Ue + - k G n Rc ++)n e - (A9) 

+ k 20 n u + k 22 n Hc + + k 21 n ne (A10) 



The equilibrium values of njj- and n H + are used as the timescales are relatively short. To conserve 
mass and charge the following local constraints are also imposed 

nu + + n H- + 2nn 2 + 2n H + = constant 
nue + n Ue + + n He ++ = constant 

n H+ - n H" + n H+ + n Hc+ + 2n He++ = n e" (AH) 
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The rate coefficients fci_ig(T) are taken from Tables A.l and A. 2 of Abel et al 199£, and the 
photoionization and photodisocciaton coefficients /C20-26 are computed from the radiation field 
(see next section). 



B. Radiative Transfer 



The photoionization and photodisocciation &20-26 coefficients are computed from the radiation 
field via 

Ui(v)J{v,r) 



ki = 47T 



hv 



-dv 



20, . . . , 26 



(Bl) 



where 020-26 (^) are taken from Table A. 3 of Abel et al 199£. The mean flux at a distance r from 
the center is computed from 



+1 



J{v,r) = iJ(^,r ha i ) / exp[-r(^,r,/x)]d//, 



(B2) 



where 



r(v,r,fjL) = ^2ai(u)Ni(r,fjL) , i = 20, ... ,26 



(B3) 



and the cross-section indices and the column densities map as follows {20, 21, 22, 23, 24, 25, 26} 
{H, He + , He, H~ , H^, H^, H2}. The column density of each species along each line of sight is 



Ni(r,n) = / rii(r (x))dx 

Jo 



(B4) 



where l(r, r^aio A 4 ) = rfi + yr^ alo — r 2 (l — /x 2 ), r'(x) = \/r 2 + x 2 — 2rxfi, and fi = cos 6 (see Figure 
9). For the impinging background UV field, we choose J{v, rh a io) = Joi^/^Ly) 01 and unless 



otherwise stated a = 1. The final coefficient for the photodissociation of H2 is taken from Draine 



fc Bertoldi 1996, which attempts to include the effects of self-shielding 



&27 (r) 

Cpump 

X 

f shieldir^) 
fshield(r, /i) 



0.15 Cpump X fshield{ r ) j 

3 x 10~ 10 sec^ , 
[J(y, r h£do )u/c] 



J c/i/=1000A 



4 x 10~ 14 erg cm -3 
1 

fshieid(r,n)dn , 



1 



j x -0.75 
10 14 cm- 2 



x < 1 
x > 1 



(B5) 
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C. Numerical Algorithm 

The following is a step-by-step outline of the code: 

1. Initialization. 

(a) Choose z v , compute Mh a i and rhaio f° r 1 _(J perturbations. 

(b) Compute M gas based on Mj eans and Meondi- 

(c) Compute pbm (f) using Burkert profile. 

(d) Choose J . Initialize T(r) ~ 20,000 °K. 

(e) Choose initial mass fractions for H, H + , H~,H 2 , H^ , ILJ, He, He + , He + +, and e~. 

2. Compute number of time steps to integrate over this value of Jo based on heating/cooling 
and chemical time scales. 

3. Compute hydrostatic equilibrium. Once the total amount of gas in the halo is specified the 
hydrostatic equilbribum gas profile can be computed by setting the gas temperature profile 
T(r). If we assume that the potential is DM dominated throughout, then 



Pgas (r)T(r) = 
p gas (0)T(0) 6XP 



p GM DM (r') 



k B T(r>) r' 2 



(CI) 



where p gas (0) is determined from the normalization M gas = J rhal ° Airr' 2 p gas (r')dr' . 



(a) Compute p gas (^)/Pgas(0) so that it is in hydrostatic equilibrium with T(r) and the DM 
potential set by /9dm(0- 

(b) Use M gas to compute p gas (0). 

(c) At each radius recompute the number densities of each species while preserving relative 
fractions. 

4. Radiative transfer. 

(a) Compute optical depth at each radius along each line of sight. 

(b) Solve for the internal radiation field J(u,r). 

(c) Compute photoionization and photodissociation rates. 

5. Non-equilibrium chemistry. 

(a) Interpolate temperature dependent rate coefficients. 

(b) Integrate non-equilibrium reaction equations. 

6. Heating and cooling. 
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(a) Compute heating/cooling due to photoionization, collisional ionization, recombination, 
collisional excitation, bremsstrahlung, Compton and H2 cooling. 

(b) Integrate and solve for T(r). 

7. If done with all the time steps for this value of Jo, decrement Jo and go to step 2. 

The result of this method is the a solution which represents the quasi-hydrostatic equilibrium 
state of the gas, which in the absence of hydrodynamic heating processes can be thought of as the 
asymptotic solution of the spherical collapse of gas into a static potential. 

At the heart of our code are the various arrays use to store the relavent data. The gas 
dynamics consists of the arrays to store the DM, gas and temperature profiles: /Odm(^), Pgas(r), 
T(r). All the radially dependent arrays are evaluated on grid rj with N r points. The chemistry 
data consist of one number density array, rii(r), for each species. 

The radiative transfer data is composed of several arrays, beginning with the array 

A jj , k = l(r j ,r , i,fik)-l(rj,r i,fj, k ) (C2) 

3 + 2 J 2 

which contains the path length spent in a shell ry of a ray traveling to a shell centered on rj 
at an angle \i k (see Figure 9). The ability to pre-compute these path lengths is critical to the 
efficient computation of the column density array for each species along each path at each radius 
N(r,[x). Using the column density and cross-section arrays ct20-26(^) the flux array can be 
computed J{y,r) and subsequently the /c2o-26- One additional complication arises from the sharp 
photoionization boundaries in the cross sections, which require very high values of N v to resolve 
accurately. However, this problem can be solved by placing extra high resolution grid points at 
these boundaries. We found 7 additional points at each boundary provided sufficient accuracy. 

The convergence properties of all the various grid parameters were tested, as well as the time 
step parameters. The values we selected were N r = 200, N„ = 20, N v = 200, which were twice the 
minimum necessary to provide accurate results. The v grid was based on a logarithmic scale with 
0.74 eV < hv < 0.74 x 10 4 eV; increasing the range of v did not alter the results. 
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Fig. 1. — Circular velocity as a function of virialization redshift for 1-a perturbations as computed 
from Press-Schecter theory using standard CDM. 
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Fig. 2. — Mhaio, ^fe-^haio, Mj eans , and Meondi for 1-a perturbations collapsing out of a 20,000 °K 
IGM. 
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Fig. 3. — Evolution of the gas density profile (top panel) and gas temperature profile (bottom panel) 
of a typical object (v c 30 kms, z v = 3, Mh a io = 2.8 x 10 9 Mq, M gas = 1.5 x 10 7 Mq) from higher 
flux (purple = 1.0 J21) to lower flux (red = 0.14 J21), where J21 = 10 -21 erg s _1 cm -2 ster -1 Hz~l. 
As the flux decreases the behavior of the core is characterized by three phases: H + , H and H2. In 
the H + phase the flux is high enough to heat the core to T core ~ 20, 000 °K. In the H phase the 
flux does not penetrate to the center and a distinct H core develops with T core ~ 10, 000 °K. In the 
H2 phase the flux is sufficiently weak that self-shielding allows H2 to form, which rapidly increases 
the cooling rate and causes the gas in the core to collapse. The transitions between phases occurs 
abruptly at Jh ~ 0.46 J21 and Jh 2 ~ 0-14 ^21 ■ 
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Fig. 4. — For each object there are critical values of the radiation flux which mark the upper 
phase transition from H + — > H (Jh) and the lower phase transition from H — > H2 («/h 2 )- The 
left panel shows the flux as a function of the central temperature for several different objects 
(z v = 4.0, 3.5, 3.0, 2.5, 2.0). We define J H to be the point where T corc drops below 12,000 °K 
(right dotted line). Likewise, Jh 2 is the point where T core drops below 8,000 °K (left dotted line). 
The left panel shows Jh,h 2 as function of the virial redshift. The straight lines are the best fits: 
ln[J H ,H 2 /J 2 i] = flH,H 2 - b u ,H 2 z v , where a H ,H 2 = (17.9, 17.3), and 6 H ,H 2 = (6.2,6.4) 
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Fig. 5. — This figure shows the collapse redshift, z v , the redshift of H formation, zr, and the 
redshift of H2 formation, function of the circular velocity of v c . Thus, a 1-a perturbation 

with v c = 30 km s _1 collapses at z ~ 3 (top arrow), but is ionized until z ~ 2.5 (middle arrow) 
and cannot form molecular hydrogen until z ~ 1.5 (bottom arrow). The ionizing field evolves as 
J(z) = Jo[(l + • z )/4] 4 , for z < 3 and J(z) = Jo> fo r z > 3 with the spectral index held fixed at 
a = 1. This figure is for a flux amplitude of Jq = 10~ . 
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Fig. 8. — Evolution of the column density profiles of H (top panel) and He + (bottom panel) of the 
same object shown in Figure 3. 
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